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CROSS-REFERENCES TO RELATED APPLICATIONS 



The present application is a continuation-in-part of United States Patent Application Ser. 
No. 09/805,422 (now United States Patent *******) filed on 13 March 2001. The 
5 present application is related to United States Patent Application Ser. No. 09/433,446 
(now United States Patent *******) filed on November 4, 1999. 

BACKGROUND OF THE INVENTION 

Field of the Invention 

[0001] The present invention relates generally to geophysical exploration and more 
10 particularly to methods for accurately estimating fluid and overburden pressures in the 
earth's subsurface on local, prospect and basin-wide scales. 

Background of the art 

[0002] During drilling of a borehole, drilling fluids, usually referred to as "mud," are 
15 circulated in the borehole to cool and lubricate the drill bit, flush cuttings from the 

bottom of the hole, carry cuttings to the surface, and balance formation pressures 
encountered by the borehole. It is desirable to keep rotary drilling mud weights as light 
as possible, but above the formation pore fluid pressure, to most economically penetrate 
the earth; heavier muds may break down rocks penetrated by the borehole and thereby 
20 cause loss of mud. Mud weight is carefully monitored and may be increased during 

drilling operations to compensate for expected increases in the formation fluid pressure. 
In some areas, however, there may be unexpected abnormal increases in pressure, with 
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depth such that mud weight does not compensate pressure; the result can be blowout of 
the well. 

[0003] Normal pressures refer to formation pressures that are approximately equal to the 
5 hydrostatic head of a column of water of equal depth. If the formation were opened to 
the atmosphere, a column of water from the ground surface to the subsurface formation 
depth could balance the formation pressure. In many sedimentary basins, shallow 
predominantly sandy formations contain fluids that are under hydrostatic pressure. 

10 [0004] At a number of offshore locations, abnormally high pore pressures have been 

found even at relatively shallow sub-sea bottom depths ( less than about 1500 meters). 
This could occur if a sand body containing large amounts of water is covered by silt or 
clay and buried. The dewatering of clays may result in the formation of relatively 
impermeable shale layers that slow down the expulsion of water from the underlying 

15 sand. The result of this is that the sand may retain high amounts of fluid and the pore 

pressure in the sand exceeds that which would normally be expected from hydrostatic 
considerations alone, i.e., the fluid pressure exceeds that which would be expected for a 
column of water of equivalent height. This phenomenon of overpressuring is well known 
to those versed in the art and is commonly referred to as "geopressure." 

20 

[0005] It is desirable to set casing in a borehole immediately above the top of a 
geopressured formation and then to increase mud weight for pressure control during 
further drilling. Setting a casing string which spans normal or low pressure formations 
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permits the use of very heavy drilling muds without risking breaking down of borehole 
walls and subsequent lost mud in the shallower interval. On the other hand, should 
substitution of heavy drilling mud be delayed until the drill bit has penetrated a 
permeable overpressurized formation (e.g., sandstone), loss of well control or blowout 
may occur. 

[0006] In areas where there is reason to suspect existence of such high pressure 
formations, various techniques have been followed in attempts to locate such geopressure 
zones. For example, acoustic or electric logs have been run repeatedly after short 
intervals of borehole have been drilled or are acquired using measurement-while-drilling 
techniques, and a plot of acoustic velocity or electrical resistance or conductivity as a 
function of depth has been made. Abnormal variations of acoustic velocity and/or 
electrical properties obtained by logging may indicate that the borehole has penetrated a 
zone of increasing formation pressure. Such techniques are very expensive and time- 
consuming and cannot predict what pressures will be encountered ahead of the bit. 

[0007] Several methods are known in the art for estimating pore pressures in formations, 
using well log data and also from seismic survey information. One such method is well 
known in the art as the "Eaton" method, and is described at Eaton, "The Equation for 
Geopressure Prediction from Well Logs" SPE 5544 (Society of Petroleum Engineers of 
AIME, 1975). The Eaton method of determining pore pressures begins with 
determination of the so-called "normal compaction trend line" based upon sonic, 
resistivity, conductivity, or d-exponent data obtained by way of well logs. The normal 



compaction trend line corresponds to the increase in formation density (indicated by 
sonic travel time, resistivity or conductivity) that would be expected as a function of 
increasing depth due to the increasing hydrostatic pressure that forces fluids out from the 
formations and thus decreases the sonic travel time (increases the velocity), assuming the 
5 absence of geopressure. This normal compaction trend line may be determined solely 

from the sonic travel time, conductivity, or resistivity well log data, or may be adjusted to 
reflect extrinsic knowledge about the particular formations of interest. The Eaton method 
calculates pore pressure by correlating the measured density information, expressed as an 
overburden gradient over depth, to deviations in measured sonic travel time, (or electrical 
10 resistivity or conductivity) from the normal compaction trend line at specific depths. The 
pore pressure calculated from the Eaton equations has been determined to be quite 
accurate, and is widely used in conventional well planning. 
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[0008] Specifically, the Eaton method determines a pressure gradient according to the 
relation 



G =G 0 -(G 0 -G„) 



-|3 



normal 



At 



observed 



(1) 



where G P is the formation pressure gradient (psi/ft), G 0 is the overburden gradient, G N is 
the normal gradient, At normal is the normal transit time and At observed is the observed transit 
time. 
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[0009] However, application of the Eaton method has been limited to the immediate 
locations of existing wells, as it depends on well log data. It is of course desirable to 
estimate pore pressure at locations at the sites of proposed new wells, and thus 
away from currently existing wells, particularly to identify locations at which production 
will be acceptable at a low drilling cost (e.g., minimal use of intermediate casing). In 
addition, knowledge of pore pressure at locations away from existing wells enables 
intelligent deviated or offset drilling, for example to avoid overpressurized zones. 

[0010] Kan (US Patent 5,130,949) teaches a method in which seismic data is combined 
with well log data to generate a two- dimensional geopressure prediction display; this 
permits deviated and horizontal well planning plus lithology detection. Shale fraction 
analysis, compaction trend, and seismic velocity may be automatically or interactively 
generated on a computer work station with graphics displays to avoid anomalous results. 
Corrections to velocity predictions by check shots or VSP, and translation of trend curves 
for laterally offset areas increases accuracy of the geopressure predictions. In particular, 
Kan '949 determines the transit time from sonic logs for compressional waves in 
predominantly shaly sections and expresses the pore pressure gradient (PPG) in terms of 
the transit time departure from the compaction trend line, 6At, as 

PPG = 0.465+ C x (8At)+C 2 (8bt) 2 m 



[0011] Coefficient C, typically varies from 0.002 to 0.02 if the transit time is expressed 
in microseconds per foot and the PPG is measured in psi/ft. C 2 may be positive or 



negative. 

[0012] Kan '949 also teaches the use of vertical seismic profile (VSP) data for 
calibration of the sonic log data. Kan (US patent 5,343,440) and Weakley (US Patent 5, 
5 128,866) further teach the use of coherency analysis of surface seismic data for 

determination of interval velocities. 

[0013] Scott (US Patent 5,081,612) teaches a variation of the Eaton method in which an 
equation of the form 

10 V c = V x (\- a x L- aj + a 3 P) (3) 

where a„ a 2 and a 3 are constants, V, is a constant, V c is calculated velocity, L is the 
lithology of the formation, cj) is the porosity and P is the effective pressure (difference 
between the overburden pressure and the formation fluid pressure). The compaction of 
the sediments is governed by an equation of the form 

15 (f> = fae"** (4) 



[0014] A reference model for the sedimentary basin is developed assuming compaction 
under hydrostatic pore pressure. A reference effective pressure and a reference velocity 
profile are obtained. An iterative procedure is used in which the lithology may be varied 
20 with depth and the reference velocity profile is compared to a velocity profile obtained 
from seismic data. 
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[0015] In addition to undercompaction, Bowers (US Patent 5,200,929) discusses a 
second cause of overpressuring. Abnormally high pressure can also be generated by 
thermal expansion of the pore fluid ("aquathermal pressuring"), hydrocarbon maturation, 
charging from other zones, and expulsion/expansion of intergranular water during clay 
diagenesis. With these mechanisms, overpressure results from the rock matrix 
constraining expansion of the pore fluid. Unlike undercompaction, fluid expansion can 
cause the pore fluid pressure to increase at a faster rate than the overburden stress. When 
this occurs, the effective stress decreases as burial continues. The formation is said to be 
"unloading." Since sonic velocity is a function of the effective stress, the velocity also 
decreases and a "velocity reversal zone" develops. A velocity reversal zone is an interval 
on a graph depicting sonic velocity as a function of depth along a well in which the 
velocities are all less than the value at some shallower depth. 

[0016] A large portion of the porosity loss that occurs during compaction is permanent; it 
remains "locked in" even when the effective stress is reduced by fluid expansion. A 
formation that has experienced a greater effective stress than its current value will be 
more compacted and have a higher velocity than a formation that has not. Consequently, 
the relationship between sonic velocity and effective stress is no longer unique when 
unloading occurs. In other words, for every effective stress, there is no longer one unique 
sonic velocity. The sonic velocity follows a different, faster velocity-effective stress 
relationship during unloading than it did when the effective stress was building. This 
lack of uniqueness is called "hysteresis." Since fluid expansion causes unloading, while 
undercompaction does not, hysteresis effects make the sonic velocity less responsive to 
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overpressure generated by fluid expansion. As a result, the pore fluid pressure 
corresponding to a given sonic velocity at given depth within a velocity reversal zone can 
be significantly greater if the overpressure was caused by fluid expansion rather than 
undercompaction. Therefore, the sonic velocity of an overpressured formation is 
indirectly dependent upon both the magnitude and the cause of overpressure. To account 
for different causes of overpressuring, Bowers teaches the use of two different velocity- 
effective stress relations: one relation applies when the current effective stress is the 
highest ever experienced by a subterranean formation and a second relation that accounts 
for hysteresis effects is used when the effective stress has been reduced. Pore fluid 
pressure is found by subtracting the computed effective stress from the overburden stress. 
Bowers uses a relationship of the form 



for the effect of unloading. In eq. (5) a, max is the maximum stress to which the rock has 
been subjected. The unloading curve parameter U is a measure of how plastic the 
sediment is, with U = 1 and U = *> defining the two limiting cases. For U = 1, the 
unloading curve is the same as the loading curve whereas for U = <*>) the velocity 
remains fixed at a value V max determined by the stress , max . 




(5) 



SUMMARY OF THE INVENTION 
[0017] The invention provides a methodology, process and computer software for the 



prediction of fluid and rock pressures in the subsurface using geophysical and geological 
data. The method includes techniques for velocity analysis from seismic data that are 
used to drive the pressure prediction, as well as an integrated approach to deriving 
pressure data that is new and novel in nature and scope. The invention addresses the 
prediction of pressure information for three scales of analysis including (1) basin-scale 
(10-500 km spatial lengths) analysis of hydrocarbon systems, (2) prospect-scale (1-10 km 
spatial length) analysis of fluid flow that can be used to analyze fluid movement in 
localized areas, and (3) prediction of pressure conditions at a specific location (0-1 km 
spatial length) where a well is to be drilled. The results of the prediction can be utilized 
in a range of other applications that address the fundamental behavior of hydrocarbon 
systems and can improve the ability to find commercial quantities of hydrocarbons in the 
subsurface. The results can also be used to design and optimize well planning. 

[0018] The technique and computer software estimates fluid pressure (from overburden 
stress and effective stress where effective stress is a function of seismic velocity), 
fracture pressure gradient, overburden pressure (from integration of density logs, using 
the Gardner-Gardner-Gregory relations with velocity estimates from sonic logs or 
seismically derived seismic velocities, or a combination of seismic data and 
Gravity/Magnetic data ), effective stress (from seismic data) and porosity (-from seismic 
data) for well planning, basin-scale pressure fields and prospect-scale pressure fields, 
and also can generate predictions and interpretations that are applicable to: 

-hydrocarbon maturation 

-fluid migration 
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-lateral and vertical seal rock integrity 
-reservoir-specific lateral pressure changes 
-fault-seal properties 

-effects of hydrocarbon reservoirs on pressure prediction 

-effects of anomalous lithologic intervals on pressure prediction and overburden 

estimation 

[0019] The software provides an operating environment in which all data pertinent to 
predicting subsurface rock and fluid pressures can be stored, displayed in tabular and 
graphical forms, analyzed, calibrated and used in the prediction process. 

[0020] The software provides the user with three pressure prediction methods, the Eaton 
method, the effective stress method, and the equivalent depth method , and allows a 
range of curve fitting options for each method including linear, power law, exponential, 
and quadratic functions. 

[0021] The computer software allows large amounts of seismic velocity data to be 
processed and analyzed with the results being displayed as a color underlayment on a 
display of a migrated or stacked seismic section. The software allows the results from an 
entire 2-D seismic line, a set of multiple 2-D seismic lines, or a 3-D seismic volume to be 
displayed with the related seismic and well data. In addition to migrated or stacked 
seismic section displays, the software allows displays based on velocity data, acoustic 
impedance data, density data, and other attributes such as frequency displays, amplitude 
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displays, phase displays, and other common derivatives created from analysis and 
inversion of the seismic data. In particular, certain displays such as acoustic impedance, 
which are generated from post-stack inversion and pre-stack inversion, can be used in the 
calibration and prediction steps of pressure prediction because of their usefulness as 
indicators of lithology variation. Likewise, amplitude and frequency displays can be 
used in the calibration and prediction steps because they identify anomalous zones where 
hydrocarbons are present in the subsurface which can cause anomalies in the data that 
would result in erroneous pressure estimates. 

[0022] The computer software allows calibration functions and prediction computations 
to be constrained by geologic boundaries such as horizons and faults so that multiple 
calibrations can be applied to areas with complex geological loading histories. 

[0023] The technique and computer software use multiple data types and handle them in 
an integrated fashion. The data include seismic and well-based velocity information 
(such as from VSP data and sonic logs), geologic boundaries such as horizons and faults 
mapped from the seismic data, pressure data from wellbores including formation fluid 
pressures and fracture pressures based on Leak Off Tests (LOT's), well log data, and 
digital seismic data in standard industry formats. LOT's are described further below/. 

[0024] The computer software has the ability to read data in generic text formats and in 
standard formats such as would be generated by well logging contractors, seismic 
processing contractors, or commercial seismic interpretation systems. The data can be 
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displayed in various tabular and graphical forms. The software provides display 
capability that includes interactive viewing and analysis of multiple data types 
simultaneously so that the user can evaluate all of these data at once. Related data can be 
calibrated using a variety of functional equations that are appropriate for pressure 
prediction. 

[0025] The technique and software contains a calibration module that includes a method 
and ability to display, edit, and datum a variety of well logs, and fit calibration curves to 
the logs or estimate equation coefficients and exponents for a specific prediction method 
(such as the Eaton method, the Bowers method, the effective depth method). The 
displays include depth versus log value displays along with cross-plot displays for 
various log parameters and plots of velocity versus effective stress. The depth and 
crossplot displays are linked interactively to a scrollable coefficient display window that 
allows the user to modify the coefficients of the equation and observe the resulting 
change in the graphical display in real time. The calibrations can be performed using all 
of the methods and equations known to experts versed in the art including but not limited 
to the Eaton method, the effective stress method and the equivalent depth method. The 
calibration module includes a method and ability for using a variety of mathematical fits 
including linear, exponential, power law and quadratic forms that encompass the full 
range of possible mathematical forms for pressure relationships in the subsurface. The 
calibration module allows the calibration to equation form of effective stress versus 
velocity, velocity versus density, velocity versus porosity, density versus porosity, 
overburden versus depth, and fracture gradient versus depth. These equations can be 
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defined for any subset of a set of well log data that the user selects, for example, on the 
basis of a lithological or other discriminator applied to the data, or by selecting a specific 
geological interval that is defined by certain characteristics and specific mapped geologic 
boundaries. The calibration of fracture gradient can be performed using LOT data and 
one of three methods which include (i) fitting a function to the LOT data from available 
local well control; (ii) determining a function based on a percentage of the overburden 
stress that honors regional data; or (iii) applying an appropriate stress ratio (KJ as 
defined by Matthews and Kelly ("How to Predict Formation Pressure and Fracture 
Gradient" (Oil and Gas Journal, v. 5, no. 8, 1967)). 

[0026] The calibration module also allows the calibration of unloading hysteresis for 
zones where the velocity has decreased due to secondary pressure conditions. This is 
determined by fitting the velocity and stress data for the unloaded zone to a relationship 
for unloading 

[0027] The calibration module includes interactive graphical displays of velocity, pore 
pressure, effective stress, overburden and fracture gradient versus depth and or time with 
horizons and faults posted on the displays. 

[0028] The technique and software contains a prediction module that includes a method 
and ability to display seismic velocity data along with displays of stacked or migrated 
seismic data with horizons, faults and an attribute of choice from a range of velocity and 
pressure data types plus the same calibration displays above. 
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[0029] The prediction module applies the calibration curves or equations to the data at 
each seismic velocity location in order to calculate such attributes as pore pressure, 
effective stress, fracture gradient, overburden, and porosity, 

5 [0030] The prediction module can store in memory the calculated attributes as functions 

of time or depth for later display and analysis. 

[0031] The prediction module displays the stacked seismic data with interpreted horizons 
and faults, and pore pressure or other attribute underlayment in color during the analysis 
1 0 that can be used interactively during the prediction of pressures. 

[0032] The prediction module includes interactive graphical displays of velocity, pore 
pressure, effective stress, overburden, density, porosity and fracture gradient versus depth 
and or time with horizons and faults posted on the displays, seismic sections and base 
1 5 maps. 

[0033] The technique and software contains a method and ability to generate digital files 
of the predicted attributes in appropriate formats for other mapping packages, and 
transfer these files to a computer medium appropriate for transfer to other computer 
20 software. Prediction results can be output in time or depth for subsequent import into 
seismic interpretation systems. 

[0034] The software contains a method and has the ability to display seismic scale 
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overlay plots of pressure and other attributes and depth plots for each location analyzed. 
The predicted attributes include, but are not limited to the following: 

-pore pressure as a function of depth or time 

-effective pressure as a function of depth or time 

-overburden as a function of depth or time 

-fracture gradient as a function of depth or time 

-porosity as a function of depth or time 

-excess pressure as a function of depth or time 

-unloading pressure as a function of depth or time 
Predicted pressure-related attributes can be displayed in units of pressure (e.g., psi) or 
pressure gradient (e.g., psi/ft or ppg). 

[0035] The technique and software includes a method and ability to use velocities from 
one or more sources including but not limited to one or more of the following: 

-stacking velocity data 

-coherency inversion velocity data 

-pre-stack inversion P-wave velocity data 

-post-stack inversion P-wave velocity data 

-pre-stack inversion S-wave velocity data 

-post-stack inversion S-wave velocity data 

-shear-wave stacking velocity data 

-tomographic P-wave velocity data 

-tomographic S-wave velocity data 

16 



-VSP velocity data 
-VSP look-ahead inversion 
-sonic logs 
-dipole sonic logs 

-mode-converted shear wave velocity data 
-combined Vp and Vs inversion 

[0036] The method includes use of all of the above velocity types for pressure prediction 
with horizon-keyed constraints. The technique also claims first use of some of the 
velocity types above without horizon-keyed constraints. These include the following: 

-pre-stack inversion P-wave velocity 

-post-stack inversion P-wave velocity 

-pre-stack inversion S-wave velocity 

-post-stack inversion S-wave velocity 

-shear-wave stacking velocity 

-tomographic P-wave velocity 

-tomographic S-wave velocity 

-VSP look-ahead inversion 

-dipole sonic logs 

-mode-converted shear wave velocity 
-combined Vp and Vs inversion 
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[0037] The technique and software includes a method and ability for picking seismic 
stacking velocities using a horizon-keyed and fault-keyed approach. This means that, 
unlike conventional seismic-velocity picking methods that use the strongest semblances 
from the seismic data, the velocities are picked to honor the formation geological and 
5 structural boundaries defined by the lithological rock units and such features as faults and 
anomalous body boundaries in 2-D or 3-D. 

[0038] The technique and computer software includes a method to isolate the velocities 
for geologic intervals that contain hydrocarbons or lithologies other than those that are 

10 appropriate to use in the prediction step. This same method can be used to isolate the 
velocities of reservoirs and seal rocks that are not in equilibrium with the encasing seal 
rocks. This method requires velocity types that have sufficient resolution to accurately 
detect the distinct velocities in these zones. In this situation, the software allows the user 
to identify the zone or zones to be excluded, and the software removes this value from 

1 5 the calculations in the prediction step. 

[0039] The software also allows the data from anomalous zones to be saved to a separate 
data file that can be analyzed separately for pressures in specific zones. 

20 [0040] The method can also be applied to pressure prediction from the results of seismic 

analysis on the geologic hazard called shallow water flows as defined in U.S. Patent 
Application Ser. No. 09/433,446 filed on November 4, 1999 having the same assignee as 
the present application and the contents of which are fully incorporated herein by 
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reference. 

[0041] The method of excluding zones from the pressure analysis can be applied for any 
velocity data types that may be developed and used in other software that are already 
5 known to those versed in the art of pressure prediction. 

[0042] The technique and software includes a method and ability to change the 
overburden and fracture gradient calculated for a given location based on one of three 
schemes. In the first scheme, the analyst identifies an anomalous density zone that may 

10 be caused by the presence of either anomalous rocks or fluids, and inputs to the software 
the anomalous density value for that zone. The software inserts this new density for the 
interval into the overburden function for the location being evaluated, and re-calculates 
the overburden function using this zone in place of the densities from the reference 
function that was used originally to determine the overburden stress versus depth. This 

15 new function with the anomalous zone is used only at the location where the anomalous 

formation was identified. By doing this at each location where the anomalous zone is 
observed, the user may build an overburden and fracture gradient profile that honors the 
actual geology in 3 dimensions rather than relying on a ID density function determined 
from a single well location that does not always represent the subsurface at other 

20 locations. 

[0043] The second scheme for changing the overburden and fracture gradient is to use 
2D or 3D density models determined independently from the inversion of vector and 
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tensor gravity and vector and tensor magnetic data and to integrate these models in 2D or 
3D to predict overburden in an area where pressures are to be predicted. This can 
include, but is not limited to the inclusion of ID well density data for calibration of the 
gravity and magnetic data. Such methods of determination of overburden are discussed, 
for example, in United States Patent Application Ser. Nos. 09/285,570, 09/399,218, and 
09/580,863. 

[0044] The third scheme for changing the overburden and fracture gradient is to use 
density values predicted in 2D or 3D seismic data using simultaneous inversion of 
multicomponent PP and PS data or the density values predicted from pre-stack inversion. 
In this case, the seismic data is calibrated with density data from well control, and the 2D 
or 3D seismically derived density values are used at each location to calculate 
overburden and fracture gradient. 

[0045] The technique and software contains a method and ability for creating and 
applying calibration functions that are controlled by and keyed to specific geologic 
intervals in the subsurface. These are usually applied to a geologic interval of specific 
age that is observed to occur across a specific area of the subsurface and has a distinct 
pore pressure calibration that is different from intervals above and below it in the 
subsurface. The method and ability allows the user to identify interpreted horizons to 
which the calibration is keyed and limits the application of the calibration to those 
horizons. 
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[0046] The method allows the user to use more than one active calibration at one time 
and the computer software splices these multiple calibrations together in the displays so 
that the user can simultaneously produce an integrated prediction for multiple zones with 
different calibrations. This multiple-zone calibration and prediction also includes the 
application of unloading parameters that can be specified for each calibration zone. 

[0047] The technique and software includes interactive help menus that guide the 
interpreter through the process and method for performing pressure prediction. The help 
menu is divided into a methods tutorial and a hands-on help manual that defines and 
specifies what each module and function in the software actually does. 

[0048] The technique and software includes a method of designing and applying a 
transformation function to data sets that are referenced to different map coordinate 
systems, including but not limited to UTM, shot point and line number, inline number 
and cross-line number, and state plane systems, such that all data can be referenced to a 
common coordinate system. Precise design and application of the transformations are 
facilitated by real-time parameter iteration made possible by concurrent examination of a 
base map display together with access to the transformation menu. 

[0049] The technique and software includes a method of organizing data in a hierarchical 
structure. At each level of the hierarchy, one and only one set of data can be active. 
Within the active data set, multiple data elements can be selected for graphical and 
computational purposes. 
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[0050] The technique and software includes a method to automatically propagate the 
effects of interactive editing of calibration functions to those graphical displays that 
depend on the edited function. This capability includes but is not limited to real time 
updating of calculated effective stress, normal trend velocity, and pore pressure values in 
5 the calibration module in response to edits to the velocity versus effective stress 

calibration function, and updating the color underlay displays in the seismic section 
display in response to edits to parameters in the prediction module. 

BRIEF DESCRIPTION OF THE DRAWINGS 

1 0 [0052] The novel features that are believed to be characteristic of the invention, both as 

to organization and methods of operation, together with the objects and advantages 
thereof, will be better understood from the following detailed description and the 
drawings wherein the invention is illustrated by way of example for the purpose of 
illustration and description only and are not intended as a definition of the limits of the 

15 invention: 

i 

Figure 1 is an illustration of a vertical section through the earth showing a well and d 
geologic horizons including faults. 

Figure 2 is an example of a flow-stream using some of the modules of the present 
invention. 

20 Figure 3 is a diagram illustrating the identification of shales at a well location. 
Figure 4a is an example of a -cross-plot of density and porosity 

Figure 4b illustrates two different examples of calibration curves relating the density and 
porosity data of Figure 4a. 
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Figure 5 shows an example of overburden stress determination: that is useful in 
calibration. 

Figure 6 shows examples of plots of velocity and porosity as a function of depth along 
with calibration curves. 

Figure 7 shows plots of velocity, effective stress and pore pressure gradient. 

Figure 8a is a plot showing a comparison of seismically derived interval velocities and 

log velocities. 

Figure 8b is a plot showing a comparison of seismically derived interval velocities and 
filtered log velocities. 

Figure 9 is a display showing seismic data and locations where seismic velocities have 
been derived. 

Figure 10 is a display of seismically derived interval velocities at two different locations 
along with a filtered velocity log from a well. 

Figure 11 shows a comparison of densities derived from seismic velocities and a density 
log. 

Figure 12 shows sonic log velocities, effective stresses in the formation and pore 
pressure data from mud weights. 

Figure 13 shows an example of a calibration of sonic velocity and effective stress using a 
Bowers curve fit. 

Figure 14 shows effective stresses derived from seismic velocities using the calibration 
of Figure 13. 

Figure 15 shows a comparison of an Eaton calibration with a Bowers calibration. 
Figure 16 shows the use of velocity spectra in obtaining seismic velocities over defined 
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seismic interval and a comparison with prior art velocity spectra. 
Figure 17 is a schematic illustration of the effects of unloading, 
b is an illustration of the equivalent depth method. 

Figures 19a-19c illustrate overpressuring that may occur in a thin sand body as a result 
5 of rapid burial. 

Figures 20a - 20b is an example of seismic data over a geologically complex prospect 
and pressure predictions made using the present invention. 

Figure 21 similar to Figure 2, is another example of a flow-stream using some of the 
modules of the present invention 
10 Figure 22 similar to Figures 2 and 21, is another example of a flow-stream using some of 

the modules of the present invention. 

Figure 23 similar to Figures 2, 21 and 22, is yet another example of a flow-stream using 
some of the modules of the present invention. 

15 DETAILED DESCRIPTION OF THE INVENTION 

[0053] Referring now to Fig. 1, an example of a vertical section of subsurface region 1 is 
shown. Indicated in the figure is a well location 25 with a well 27 penetrating the 
subsurface. A number of faults 10, 11, 12, 13, 14, 15, 16, 17, 18 are indicated in the 
figure as well as a number of horizons 21 that correspond to geologic intervals of 

20 interest. 

[0054] In the well 27, a plurality of measurements may be made of the properties of the 
subsurface formations penetrated by the well. These typically include sonic logs that 
measure the velocity of compressional and shear velocities, density logs, gamma ray logs 

24 



that are indicative of the shale content of the formation, and resistivity logs of various 
types that measure the formation resistivity. 

[0055] In addition to these logs, a record is kept of the mud weight that is used for the 
5 drilling of the wellbore: as noted above, the mud weight is usually selected to maintain a 

slightly overbalanced condition wherein the borehole fluid pressure is slightly in excess 
of the formation fluid pressure. After the well is drilled, pressure drawdown and pressure 
buildup tests may be run at selected depths to make an accurate determination of the 
formation fluid pressure. These are indicative of the fluid pressures encountered in the 
10 well. 

[0056] A line of 2-D seismic data (not shown) may be acquired over the region. 
Alternatively, a plurality of 2-D seismic lines or a set of 3-D seismic data may be 
acquired over the surface of the earth that includes the section shown in Fig. 1 . 

15 

[0057] In one embodiment of the present invention, using the log information and/or 
seismic information in the vicinity of the well, a determination of the effective pressure 
in the formations and the overburden pressure on the formations is made. In another 
embodiment of the invention, the seismic data are used to obtain estimates of the 
20 effective stresses in the subsurface at locations away from the well. This may be done in 
the immediate vicinity of the well; on a prospect scale location (e.g. between the faults 11 
and 12); or on a basin wide scale (e.g., extrapolating to another prospect such as the 
region between the faults 13 and 18). 
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[0058] In the present invention, the interpreter has a plurality of modules that perform 
different functions to enable this process of extrapolation away from the well. 
Depending upon the nature of the data available, different combinations of these modules 
may be used by the interpreter. Rather than providing an exhaustive list of the various 
combinations that may be used by the interpreter, the description herein first provides 
examples of the processes performed by the various modules and then gives examples of 
sequences that may be used. 

[0059] Turning now to Fig, 2, an overview of one method of using the present invention 
is illustrated. The input data 51 includes logs, pressure measurements at the well, 
velocity functions characterizing the propagation of seismic waves in the subsurface, 
seismic data, seismic horizons picked on the seismic data and surface elevations (for land 
seismic data) and water depth (for marine seismic data). These are described in more 
detail below. 

[0060] The logs that may be used in the present invention include conventional sonic 
logs that measure the compressional velocity of elastic waves in the formation. These 
would be known to those versed in the art and are not described further. In addition to 
conventional sonic logs, an embodiment of the present invention also uses measurements 
of shear velocities in the formation. These may be direct measurements of shear 
velocities using devices such as a dipole logging tool. U.S. Patent 4,606,014 to Winbow 
et ah gives an example of a dipole sonic logging tool. Other examples of dipole logging 
tools are described in U.S. Patents 4,703,460 to Kukjian et al., U.S. Patent 4,862,991 to 
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Hoyle et al. and are not discussed further. Alternatively, direct measurements of shear 
velocities may also be obtained using a quadrupole or octupole logging device. Such 
devices are described in U.S. Patent 4,855,963 to Winbow et al. and U.S. Patent 
4,951,267 to Chang et al. and are not discussed further. 

[0061 J Shear velocity measurements of the formation may also be obtained by indirect 
methods. For example, U.S. Patent 4,774,693 to Winbow et al teaches a method for 
determination of shear velocities using guided waves produced in a borehole by a logging 
tool. In the present invention, shear velocities obtained by any direct or indirect method 
may be used. 

[0062] Shear velocities are an important indicator of the effective stress of a fluid filled 
rock. U.S. Patent Application Ser. No. 09/433,446 of Huffman having the same assignee 
as the present application and the contents of which are fully incorporated herein by 
reference teaches use of the relation between effective stress and the shear velocity of 
clastic sediments. An optional embodiment of the present invention uses this relationship 
to obtain an estimate of effective stress from shear velocity measurements. 

[0063] Density logs may also be used in the present invention. They may be used to 
obtain an estimate of the compressional velocity using the Gardner-Gardner-Gregory 
(GGG) relation 
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p = 023V 25 (6) 

where p-' is the density in gm/cm 3 and V is the compressional wave velocity in ft/s. An 
optional embodiment of the present invention uses a general power law relationship in 
which the exponent and the constant of proportionality may be determined from a 
combination of density logs and sonic logs. The density logs may also be used in the 
determination of the overburden stress as described below. 



[0064] Logs indicative of the shale content of the formations (such as gamma ray logs) 
may also be part of the input to the present invention. As described below, portions of 
the sonic logs corresponding to shaly intervals may be used as part of the calibration 
process. In a preferred embodiment of the invention, the shale-prone intervals are 
selected preferentially and are used to define a shale compaction trend in the calibration 
process. 

[0065] Pressure measurements may also be recorded and used including Mud Weight 
data, Repeat Formation Tester (RFT) and Modular Dynamic Tester (MDT) for estimates 
of the formation fluid pressure. The MDT is the new Schlumberger tool that does the 
same thing as the RFT tool, but also has several new features including (1) a sensor in the 
device, that tests fluid resistivity to make sure that the formation fluid being gathered is 
pristine formation fluid, (2) a larger set of sample chambers for gathering fluid, and (3) a 
pump that actually draws fluid out of the formation rather than relying on ambient flow 
to fill the chambers. Information on the fracture pressure of the formation, which is 
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determined through the use of LOT's, and is also expressed in gradient form as the 
Fracture Gradient, can also be used in the software. The LOT is a standard test that is 
well known to those versed in the art, and is used to determine the pressure at which the 
formation will begin to fail through tensile fracture. The test is commonly performed by 
increasing the fluid pressure in the wellbore after a new casing string has been put in 
place and cemented so that only a small interval of the well is tested. This constraint 
assures that the formation failure zone is clearly known from the test. As the pressure 
increases during the pumping phase of the test, it pushes on the exposed rock formations 
in the open portion at the bottom of the wellbore. As the mud pressure exceeds the 
fracture strength of the formation, fluid begins to "leak off into the formation as the 
fracture opens. The "leak off point" as it is called, is the pressure at which the 
compression curve begins to change its slope due to the loss of fluid into the formation. 
A standard LOT is performed by pumping the fluid pressures up to the leak off point, 
then shutting in the well and allowing the fluid pressures to slowly decline until the 
fractures push the fluid back into the wellbore and the fractures close. The leak off point 
has been determined to be close to the fracture initiation pressure, and this pressure value 
is commonly used for fracture pressure and fracture gradient estimation. 

[0066] The use of velocity information from seismic data has been described in prior art. 
However, the use of such velocity information in prior art is limited to conventional 
coherency spectra of gathers of common-midpoint (CMP) data. In the present invention, 
a novel aspect is the use of such conventional coherency spectra in conjunction with 
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picked horizons of the seismic data. This is illustrated below and discussed with 
reference to Fig. 16, In addition, alternative embodiments of the present invention use 
other methods for determination of velocities of seismic waves in the subsurface. 
Specifically, in land seismic prospecting, shear seismic data may be used and velocity 
spectra derived therefrom. In addition, in marine seismic prospecting, a combination of 
compressional and converted wave seismic data may be used for determination of 
compressional and shear velocities. The use of converted wave prospecting for velocity 
determination is discussed, for example, in US patent 4,881,209 to Bloomquist et al., US 
Patent 6,128,580 to Thomsen and is not discussed further here. 

[0067] Alternative embodiments of the present invention use velocities derived from 
prestack inversion of compressional wave data , post-stack inversion of compressional 
wave data, prestack and post-stack inversion of shear seismic data, joint inversion of 
compressional and converted wave seismic data and joint inversion of compressional and 
shear seismic data. An example of prestack inversion of compressional wave seismic 
data is given in United States Patent 5,583,825 to Carrazonne et al. the contents of 
which are incorporated herein by reference. Other examples of velocity determination 
using inversion techniques are given in Savic et al, Mallick et al., Connolly, and Duffaut 
et al. 

[0068] Post-stack inversion is one alternative to conventional velocity analysis that 
provides higher resolution by inverting for impedance from the reflection strength. Post- 
stack inversion allows the analyst to separate the seismic wavelet from the reflection 
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series represented by the geologic formations, and results in an estimate of residual 
impedance for each layer. Post-stack inversion can be applied using only the stacked 
seismic data, or can be calibrated with well logs, check shot surveys, VSP data and 
seismic velocity data. When calibrated properly with low frequency velocity field data, 
the analysis can be used to generate an estimate of the absolute impedance, or its 
components of velocity and density. This requires a good set of density-velocity 
relationships for the lithologies encountered in the wells, so that these two components 
can be separated effectively. This exercise is not trivial, however, because the post-stack 
inversion technique ignores the fact that offset-dependent behavior (AVO) is buried in 
the stacked response and can cause significant perturbation of the results. One way to 
overcome this limitation and also boost resolution of the results at the same time is to use 
only the near offset traces for the analysis. This is a good method to use as it provides 
higher resolution results due to the removal of far-offset data that are degraded by NMO 
stretch. Near-stack inversion also gives the most robust calibration to well logs that are 
essentially measuring the same vertical-incident information. 

[0069] In an optional embodiment of the present invention, pre-stack inversion of 
seismic data provides a higher level of accuracy and allows prediction of pressure at a 
scale that was not achievable in the past. Pre-stack inversion can estimate the P wave 
velocity, shear wave velocity, and density simultaneously by using the near-offset 
reflectivity and amplitude versus offset behavior of each reflection in the subsurface. 
This allows the interpreter to estimate the overburden and effective stress from the same 
data set. The resolution limits for pre-stack inversion are approximately at the tuning 
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thickness of the individual formations, so that pressure data can be generated for layers 
on the order of 100-200 feet thick at moderate depths in clastic basins. The biggest 
drawback to pre-stack inversion at this time other than cost is its extreme sensitivity to 
data quality. A robust inversion requires data that are relatively clean, uncontaminated 
by multiples and noise, and in areas of little structural complexity. In addition, the 
resolution and accuracy of pre-stack inversion is only as good as the quality of the 
reflection events in the data. If there are pressure cells that do not have reflections 
associated with them, no velocity or reflection-amplitude technique including inversion 
will identify those zones. 

[0070] The use of calibrated pre-stack inversion, especially where multicomponent data 
are available, allows the estimation of density and velocity simultaneously along a line. 
This makes possible predictions that take into account the lateral variations in density 
that accompany changes in pressure. At present, most methods of pressure prediction 
assume that the overburden from the control well represents overburden globally, which 
is false in general. Pre-stack inversion helps to remedy this problem by predicting the 
density from the seismic data and allowing the pressure interpreter to make judgements 
about the density using the seismic and well data and his own intuition. 

[0071] Pre-stack inversion makes possible isolation of velocities for individual sand 
packages so that the interpreter can determine where disequilibrium may exist between 
the sand-bearing formations and massive shales, and to isolate the velocity and density 

32 



effect of hydrocarbon-bearing reservoirs on the velocity field around them. At present, 
most methods lump these effects into thicker stratigraphic intervals that have a single 
velocity attached to them and hide the effect. These errors can cause predictions to 
overestimate or underestimate pressures significantly which leads to less effective well 
planning. 

[0072] Velocity information for the present invention may also be obtained using cross- 
well tomography. This4s generally has poorer resolution in a horizontal direction but 
may have higher resolution in a vertical direction than analysis of velocity spectra of 
surface seismic data. Depending upon the type of sources and detectors used, 
compressional and/or shear velocities may be obtained. Such tomographic methods 
would be known to those versed in the art and are not discussed further. 

[0073] Velocity information in a limited region in the vicinity of a wellbore may also be 
obtained using Vertical Seismic Profiling (VSP). These may be for P- wave velocities or 
for S-wave velocities and may be either in a conventional VSP or in a look-ahead mode 
wherein velocities are obtained pertaining to formations ahead of a drilled well. Such 
methods would be known to those versed in the art and are not discussed further. 

[0074] Returning now to Fig. 2, one of the steps that may be used in calibration is to 
identify the shale formations in the well 53. The reason for this is that clay and shale 
tend to show the effects of compaction better than sand. Fig. 3 shows an example of the 
identification of the shales in the subsurface. The left panel 101 shows a sonic log 111, 
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the center panel 103 shows an example of a porosity log 113 and the right panel 105 
shows an example of a gamma ray log 115, The porosity log 113 is obtained external to 
the present invention from the density log. The density porosity is preferable to a 
neutron porosity. The gamma ray log 115 is a good indicator of shaliness in clastic 
5 sediments because shales have a higher potassium content than sands as a result of which 

shales have a higher gamma ray count. Indicated in the right panel 105 is a shale line 
117 that may be determined interactively by the user of the invention. In a preferred 
embodiment of the invention, the shale line may be represented by a piecewise linear 
function. Only velocity and porosity samples that correspond to depths where the gamma 
10 ray count is to the right of the shale line 117 are used in the calibration process. An 

optional plot in the present invention (not shown) provides the interpreter with a color 
display of the selected and deselected portions of the various logs. 

[0075] A useful aspect of the present invention is the ability to cross-plot different 
15 properties. An example of this is shown in Fig, 4a that is a cross-plot of density along 

the abscissa 131 and porosity as the ordinate 133. The present invention provides the 
user with the ability to interactively edit the data by drawing a polygon such as 135 
enabling the selection of points within the polygon and deselecting points outside the 
polygon. As with other displays, the selected and deselected points appear in different 
20 colors. 

[0076] Fig. 4b shows the data of Fig. 4a with two different curves 137 and 139 
superposed thereon. In the example of fig. 4b, the curve 137 corresponds to the GGG 

34 



relation whereas the curve 139 is a different fit determined interactively by the 
interpreter. The interpreter has a variety of fits available, including power laws, 
exponentials and polynomials. One of the useful aspects of the present invention is that 
in the interactive displays, an active curve or data points appear in one color (e.g., green) 
5 and the results of user interaction appear in a second color (e.g., red). This enables the 

interpreter to easily compare the current iteration with what may have been determined 
earlier. Different curves may be derived in different regions of the plot. 

[0077] Cross-plots and calibrations such as shown in Figs. 4a - 4b are also indicated at 
10 59 in Fig. 2: different calibrations may be derived by the interpreter for the density- 
velocity relation and the porosity-velocity relation. 

[0078] Still referring to Fig, 2, the present invention provides the interpreter with the 
ability to determine the overburden stress 55. This is illustrated in Fig. 5. Two panels of 
15 data are shown. The right panel 151 is a plot of density (abscissa) vs. the depth 161 

(ordinate). The density values from the density log are shown 162: these may or may not 
have been preselected as described above. Also shown is an example of a power law 
curve 163 that starts with a density of 1.4 at the ocean bottom. Additional points 164 
may be added by the interpreter. 

20 



[0079] Shown in the right panel 153 is the determined overburden stress. Two curves are 
shown: 167 corresponds to the curve of density-depth values between the sea floor and 
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the top of the log shown in the left panel 151 while the curve 165 corresponds to a 
uniform density of 2.0 from the sea floor to the top of the log. As in other displays, the 
interpreter has the capability of changing the parameters of the curve fit and immediately 
seeing the effect this has on the determined overburden stress. 

[0080] Another display that is important is shown in Fig. 6. The right panel 203 is a 
cross-plot of density against velocity. This is the same data that is in Fig. 4b; however, 
two types of symbols are used in the plot. The light colored symbols are for the shale 
samples selected in Fig. 3 while the dark symbols are used to plot the remaining values 
from the logs of Fig. 3. The left panel 201 is a cross-plot of velocity against porosity, 
with all the data points shown. Also shown in the right panel 203 are two different curve 
fits to the data: the curve 205 corresponds to a density constrained to be 1.4 at the ocean 
bottom (velocity =5000 ft/s) while the curve 207 is an unconstrained fit. 

[0081] Turning now to Fig. 7, another display that is useful in performing the calibration 
is shown. The left panel 225 is a plot of the velocity log 231. In the actual screen 
display, the shales and the non-shales would appear in different colors. The right panel 
229 is a plot of the overburden stress 233 on a pore-pressure gradient (PPG) scale. Also 
displayed in the right panel are mud weight data from the scout data and the mud log 235, 
237. The center panel 227 on a scale of stress shows the overburden stress 239 and the 
effective stress 241. The effective stress in this example is based on an assumption that 
the fluid pressure is hydrostatic. The difference in the appearance of the overburden 
stress in the center 227 and the right panel 229 is because 227 is in pressure and 229 is in 
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pressure gradient units. The reason the effective stress appears as a fat line is that it 
comprises diamond shaped plot symbols every 50 ft. 

[0082] Those versed in the art would recognize that sonic log velocities do not 
5 necessarily agree with seismically determined velocities in the immediate vicinity of the 

well. This is brought out in Figs. 8a and 8b. Figure 8a is a plot of the seismically 
derived interval velocity 255 and the sonic log velocity 257. The abscissa 251 is the 
velocity and the ordinate 253 is depth. Fig. 8b is similar to Fig. 8a except that the .sonic 
log displayed is a filtered log 257'. They clearly show that the seismic velocities are 
1 0 higher than the sonic logs. 

[0083] The reason for this discrepancy becomes clear upon examination of Fig. 9 which 
is a display of the seismic section with the well location 271 indicated, along with the 
location 273 where the velocity used in the display of Figs. 8a, 8b was derived. The well 

1 5 is seen to be on a flank of a salt dome 279 with steeply dipping sediments on its flanks. 

In addition to the dips, the shallow seismic data at location 273 is noisy, possibly due to 
gas leakage from a reservoir. Fig. 10 shows a display of the filtered sonic log 257' along 
with the seismic velocity 255 at location 273 and the seismic velocity 285 from a location 
275 in Fig. 9 away from the salt dome. The velocity function 285 is in good agreement 

20 with the filtered sonic log 257\ 

[0084] Another useful display in the calibration process is a comparison of densities 
derived from seismic velocities with the density log at the well. This is shown in Fig. 11 
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where the left panel 301 is identical to Fig. 10 and the right panel shows the density log 
305 at the well, and density functions 307 and 309 derived from the seismic velocities at 
locations 273 and 271 in Fig. 9 respectively. The density function is derived from a 
density-velocity calibration such as curves 205 or 207 in Fig. 6. A similar plot may be 
made (not shown) of seismically derived porosities and well porosities. 

[0085] In the present invention, a useful check of the calibration process is to compare 
seismically derived velocities, densities and porosities with log data from another well if 
such a well is available. Such a comparison serves as an indication of how far from a 
calibration well the seismic data may be used for prediction of formation properties. If 
such a comparison is good over certain geologic intervals but not over other geologic 
intervals, this enables the interpreter to use a different calibration for different regions of 
the subsurface region 1 in Fig. 1. 

[0086] Fig. 12 is a plot similar to Fig. 7. The left panel 325 shows the sonic log 331 as a 
function of depth. The right panel 329 is a display on a pressure gradient (mud weight) 
scale of the overburden pressure gradient 333 and the mud weight observations at 
specific depths 337. The center panel is a display on a stress scale of the overburden 
stress 335 and the effective stress 337' corresponding to the observed mud weights 337. 
Using the data from the middle panel 327 a display of the effective stress calibration may 
be obtained as shown in Fig. 13. Shown in Fig. 13 is a plot of the sonic velocity as a 
function of the effective stress. The velocity data points corresponding to the mud 
weight observations 337 ' are plotted as the points 337". The horizontal scale is a stress 
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scale and a fit curve 351 is indicated on the display. As indicated in the figure, the 
interpreter has control over the type of curve fit (linear, quadratic, exponential or 
Bowers) 353 and the parameters of the fit are shown on the screen display, 

5 [0087] Using the results of the calibration curve 351 in Fig. 13, seismically derived 

velocities such as those shown in panel 301 of Fig. 11 maybe used to estimate the 
effective stress. This is shown in Fig. 14 where the left panel 352 shows the seismically 
derived velocities 357 near the well location and the calibration curve for normal 
pressure, the center panel 353 shows the estimated effective stress curve 359 along with a 

10 display of the overburden stress calibration 361. The right panel 355 shows, on a pore 

pressure gradient scale, the estimated pore pressure curve 359' along with the overburden 
stress gradient 363. As can be seen in Fig. 13, the effective stress is well below what 
would be expected under hydrostatic conditions, indicating that the entire interval is 
overpressured. A display similar to Fig 14 but not shown here may be used to shown a 

1 5 comparison between the effective stress expected on the basis of the sonic logs with the 
expected stress under normally-pressured conditions: such a comparison may be used 
prior to the prediction of effective stresses from seismically derived velocities shown in 
Fig. 14. 

20 [0088] In a manner similar to that described immediately above for the Bowers 

calibration, the present invention also provides a capability for doing a calibration based 
on the Eaton relationship. This is not described here, but the end result of an Eaton 
calibration is shown in Fig. 15 as a comparison with the Bowers calibration. As in Fig. 
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14, the left panel 375 shows a velocity plot, the middle panel 377 is a plot of effective 
stress and the right panel 379 is a pore pressure gradient plot. The middle panel shows 
the effective stress 381 from the Eaton calibration and the right panel shows the pore 
pressure gradient 381 ' from the Eaton calibration. It may be seen that the Eaton 
calibration does not match the data points for the mud weight below about 5000 ft. in 
both the panels. 

[0089] The above discussion dealt with the use of data from a well corresponding to a 
specific calibration location within the areal extent of the seismic data with or without 
seismic velocity data for the specific well location. As an alternative to deriving the 
calibration curves, the present invention includes the capability of using a previously 
determined calibration curve in the analysis of the seismic data using regional 
information or information from a location that is outside the areal extent of the seismic 
data. 

[0090] The above discussion also dealt with the use of density data from well logs to 
establish an overburden stress. United States Patent Applications Ser. Nos. 09/285,570 
filed on April 12, 1999 (now US Patent 6,278,948); 09/399,218 filed on September 17, 
1999 (now US Patent 6,424,918); and 09/580,863 filed on May 30, 2000, all having the 
same assignee as the present application and the contents of which are fully incorporated 
herein by reference, teach the use of potential fields data, including scalar, vector or 
tensor gravity and/or magnetic data, in the determination of densities of subsurface 
formations in conjunction with seismic data. The teachings of these applications may be 
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used in the present invention for obtaining densities of subsurface formations; as 
discussed in the aforesaid applications, the density model may be 1-D, 2-D, 2.5-D or 3-D. 
In addition, United States Patent Application Ser. No. 09/405,850 filed on September 24, 
1999 (now US Patent *******) and having the same assignee as the present application, 
teaches the use of potential fields data in combination with seismic data for obtaining 
estimates of overburden stresses and effective stresses in subsurface formations. The 
teachings of the '850 application may also be used in the present invention. 

[0091] As mentioned above, seismic velocities may be derived in the present invention 
using many methods. The discussion above relating to Figs. 8-15 could use velocities 
from any of these methods. In one embodiment of the invention, the velocities may be 
derived without reference to seismic horizons defined by the interpreter on a seismic 
section. In another embodiment of the invention, the seismic velocities are defined with 
reference to seismic horizons defined by the interpreter. 

[0092] One such method that has been used in the past is the use of velocity spectra for 
the determination of interval velocities of compressional waves. However, these have 
been done in prior art without reference to defined seismic horizons. Fig. 16 shows how 
velocity spectra may be obtained with referencing to defined seismic horizons. The left 
panel 401 shows conventional velocity spectra for compressional wave data. The vertical 
axis is seismic travel time and the horizontal axis is the moveout velocity. Plotted in the 
panel 401 are coherency values of the seismic data at different moveout velocities and 
times: such spectra would be well known to those versed in the art. In prior art methods, 
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peaks of the spectra such as at times indicated by 405a, 405b, 405c, . . . 405n are picked 
and seismic interval velocities derived from these peaks using known methods. As 
would be known to those versed in the art, for causes beyond the scope of the present 
invention, these peaks do not necessarily coincide with geologically meaningful seismic 
5 horizons. The right panel 403 shows the same velocity data with seismic times indicated 

as 407a, 407b, 407c 407 k that correspond to seismic horizons on a seismic section 

(not shown). In the present invention, coherency values at these times 407a, 407b ... are 
picked to define the moveout velocity as a function of depth. These do not necessarily 
correspond to local peaks of the velocity spectra. The use of horizon keyed velocity 
1 0 spectra of compressional waves is a novel aspect of the present invention. The other 

methods of determination of seismic velocities described above are believed to be novel 
with or without the use of horizon keyed picking. 

[0093] The present invention also has the capability of accounting for the effect of 
1 5 unloading of the effective stress in the subsurface formations. The effect of unloading is 

best understood with reference to Fig, 17. 

[0094] This diagram is a 3-dimensional view that demonstrates the interplay between 
velocity 503, porosity 501 and effective stress 507. The loading path starts at an 
20 effective stress of zero, and the velocity increases and porosity decreases until the 

material changes over from a Wood's Equation material to a frame-bearing clastic rock 
that can support an effective stress on the grains. The Wood's Equation portion of the 
loading path 511 occurs as the material is initially deposited and compacted near the 
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surface. Once the critical porosity such as 512 is reached, the material follows the 
primary compaction curve 513, achieving either a compacted or undercompacted state. 
If allowed to compact normally with fluid draining out of the pore spaces, a rock will 
continue up the normal loading path and velocity will increase and porosity will 
decrease. Both of these properties are dependent on the effective stress on the grains that 
are bearing the external load. If at some point the fluid is prevented from escaping, the 
rate of ascent up the normal pressure curve will decrease so that the rock has a lower 
velocity and effective stress than would be expected at normal pressure conditions at a 
given depth of burial. This condition is known as undercompaction or compaction 
disequilibrium. The key to understanding undercompaction is to recognize that a rock 
under these conditions still remains on the normal compaction trend, only it is not as 
compacted as you would expect it to be at that depth of burial under normal hydrostatic 
pressure. 

[0095] Unlike undercompaction, a rock subjected to secondary pressure (also called 
unloading) cannot stay on the normal compaction curve. When fluid is pumped into a 
rock or expands within the pore spaces in the rock, the compaction process is arrested 
and the rock begins to display a form of hysteresis behavior in velocity-effective stress 
space. When this occurs, the porosity essentially does not change except for some minor 
elastic rebound (Moos and Zwart, 1998), and the velocity behavior is strictly controlled 
by the contact area and the grain-to-grain contact stresses in the rock. Because there is 
essentially no porosity change, the net effect is to flatten out the velocity-effective stress 
trend and produce an unloading trend that is different from the compaction trend. The 
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unloading curve must start from the velocity-porosity-effective stress point on the 
primary compaction curve where the unloading begins. This is why unloading 515 
always starts from a porosity-velocity-effective stress point on the primary compaction 
curve. Note that the unloading paths occur essentially in the velocity-effective stress 
5 plane as the porosity decrease associated with compaction is arrested during unloading 

and very little elastic rebound (less than 1 porosity unit) occurs during the unloading 
process. As the effective stress decreases due to higher fluid pressures at fixed 
overburden, the velocity decreases in direct relation to the stress change. Once a rock is 
on an unloading path, the rock doesn't change porosity unless other phenomena such as 

10 diagenesis or cementation are occurring concurrently with the pressure changes. For the 

rock to begin compacting again, the secondary pressures must first bleed off, or the 
overburden must increase sufficiently by additional sediment loading to counterbalance 
the secondary fluid pressures that were added to the rockmass and increase the effective 
stress. In either case, the rock will respond to the change in effective stress and will 

15 move back up the unloading path 515 until it contacts the normal compaction curve 

again. Once the effective stress has exceeded the value where unloading began, the rock, 
can begin to compact again. If the stress never reaches this level, the rock will remain on 
the unloading path indefinitely. It is important to recognize in this context that the 
normal compaction curve is also the maximum compaction, maximum velocity and 

20 minimum porosity that a material can achieve at normal pressure for a given effective 
stress. 

[0096] To properly predict pressure ahead of the bit, it is essential to know not only the 
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normal compaction trend, but also the slope of the secondary pressure curve and the 
maximum stress- velocity state that was achieved before unloading began. It is important 
to recognize that the presence of two possible pressure mechanisms and a range of 
possible maximum velocities for unloading leads to a range of possible predicted 
pressures according to which mechanism is assumed to be at work and where unloading 
began. For any velocity, there are a range of possible pressures that are a function of the 
normal trend, the maximum velocity attained by the rock, and the unloading curve slope. 
In practical terms, for any observed velocity value, the minimum pressure case is 
represented by the normal trend curve (equivalent depth of burial) and the maximum 
pressure case is represented by the greatest reasonable maximum velocity on the normal 
trend and the slope of the unloading curve from that point back to the observed velocity 
value. Thus, it is imperative that the pressure prediction expert be aware of both causes 
of pressure and also recognize when and how to apply unloading. 

[0097] Turning now to Fig. 18, the equivalent depth method is illustrated. This method 
is applicable only for undercompaction and not for modeling the effect of unloading. 
Shown is a shale parameter (e.g., effective stress, velocity, etc.) plotted as a function of 
depth. A number of data points 527a, 527b . . .527i . . . 527n are shown along with the 
normal compaction trend curve 525. Due to rapid burial, the points below 527 are 
undercompacted, as indicated by the reversal of the trend data. Data from the depth 531 
have the same properties as data along the normal compaction curve at a depth 529. The 
depth 529 may be referred to as the equivalent depth to 531. 
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[0098] An isolated sand layer within a thick shale that is subjected to rapid burial may 
have very unusual stress configurations. This is illustrated in Fig.s 19a - 19c (after 
Stump et al). Consider a sand body 551 as shown in Fig. 19a that is initially in a 
horizontal position and then due to rapid burial at the right end, assumes the 
configuration shown by 551 f in Fig. 19b. Consider now the relative pressures between 
the sand and the shale at the shallow end (points 555, 553) and the deep end (points 
556,554). Normal hydrostatic and lithostatic stress distributions are indicated in Fig. 19c 
by 571 and 573 respectively. The shale 553 at the shallow end is essentially at 
hydrostatic pressure given by the point 553 f while the shale at the deep end 554 is at an 
abnormally high pressure denoted by the point 554'. (If the subsidence is rapid enough, 
the shale follows a stress line 575 parallel to the lithostatic line 573). The sand at the 
deep end will now be at a pressure denoted by 556' but due to the good permeability of 
the sand, the pressure gradient within the sand will be substantially hydrostatic and the 
shallow end of the sand will now be at a pressure denoted by 555'. As a result of this, the 
stress in the sand is greater than the stress in the adjoining shale and, if the difference is 
large enough, this can lead to a breakdown of any possible sealing strength of the sand-- 
shale interface and any hydrocarbons that may be present in the sand will leak out. 
Seismic inversion techniques discussed above (prestack, post-stack, PP, PS) are 
particularly useful in identifying trap integrity problems involving thin sand layers. The 
use of these methods to isolate the velocity of a high-permeability reservoir zone and 
then predict pressures as a function of depth or structural position is a preferred 
embodiment of the present invention. 
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[0099] The stresses determined by using the method of the present invention can thus be 
used as an indicator of trap integrity in hydrocarbon reservoirs. 

[0100] Referring now to Figs. 20a and 20b, the seismic data corresponding to Fig. 1 are 
5 shown in Fig. 20a. The position of the well 25 of Fig. 1 is shown at 625 along with the 

location of a second well 627, and the locations of the various faults and horizons shown 
in Fig. 1. Fig. 20b is a display of the interpreted fluid pressure in the region 
encompassed by the seismic data. Based on the display, fault 610 shows a major fault 
seal failure near the horizontal position indicated by 650 while the fault 612 shows a 
10 major fault seal failure near the horizontal position indicated by 652: the pressure is not 

discontinuous across the faults at the indicated location. In contrast, the figure shows a 
pressure discontinuity across the fault 613 indicative of fault seal integrity. There also 
appears to be a somewhat less serious fault leak across the fault 618. 

15 [0101] Those versed in the art would recognize that trap integrity may also be used in the 

vertical direction: a pressure discontinuity across a known reservoir caprock is indicative 
of vertical integrity of a seal. On the other hand, if pressure appears to be continuous 
across a caprock, then the trap integrity is questionable. 

20 [0102] Use of such displays on a basin wide scale is clearly useful in basin modeling. In 

addition, using geologic information, it may be possible to identify source rock intervals 
on a prospect scale or a basin-wide scale. Abnormally high fluid pressures in such source 
rock intervals are indicative of a secondary buildup of pressure caused by source rock 
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maturation and a subsequent expulsion of hydrocarbons into the rock matrix. 

[0103] The pressure displays along with knowledge of lithologies associated with 
different subsurface formations can provide information about the migration of fluids in 
the subsurface. This is an important diagnostic in prospect evaluation as a trap with a 
large "drainage" area is likely to contain more hydrocarbons than one with a small 
drainage area. The pressure displays may also be used to estimate lateral pressure 
changes within a reservoir. The pressure displays may also be used in the planning of 
drilling of wellbores: abnormally high fluid pressures indicated on seismic data would 
indicate the necessity of using higher mud weight in the drilling of a well. This may be 
used in the analysis of wellbore stability. Abnormally high fluid pressures indicated on 
seismic data near the ocean bottom are a warning of shallow intervals where blowouts 
may occur. 

[0104] Returning now to Fig. 2, the present invention also has the capability of 
producing maps 59. Any of a plurality of parameters may be displayed on a base map in 
a manner suitable for helping an interpreter analyze a prospect. For example, displays 
may be made of a parameter between any two seismic horizons, at any depth or at any 
seismic reflection time. To facilitate this display, location transforms may be applied 57 
on the well positions, the positions of seismic lines, seismic horizons, velocity functions 
etc. Methods for applying these transforms and producing a suitable map would be 
known to those versed in the art and are not discussed further. 
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[0105] Given the variety of modules and options in the present invention, by now those 
skilled in the art would have recognized that other flow charts besides the one shown in 
Fig. 2 may be used. Fig. 21 shows a flow chart of operations that are particularly 
suitable for processing large amounts of data with a minimum of human interaction. 
Data are input 701 relating to velocity functions, horizons, elevations (for land seismic 
data) and bathymetry (for marine data), and locations of seismic lines. Predetermined 
calibrations are input 703 describing the overburden stress as a function of depth and a 
relation between seismic velocity and effective stress using, e.g, the Bowers, Eaton or 
other suitable method. Optionally,-a density-velocity and porosity-velocity relations may 
be input 705. Using the velocity functions from 701 and the predetermined calibrations, 
predictions are made of parameters that may include the pore pressure, density and 
porosity using any of the methods described above. The predicted parameter values 709 
may be exported for use elsewhere 717, and/or displayed on a seismic section 715. In 
addition, maps may be produced 713 using the location transform information 711. 

[0106] Fig. 22 is a flow chart that may be used in using the basic Bowers or Eaton 
method. The input data 801 comprise logs (e.g., sonic, density, gamma ray, porosity), 
pressure measurements in a well (RFT, MDT, mud weight, LOT), velocity functions, 
seismic horizons, elevations (for land data) and bathymetry (for marine data) and seismic 
data. Shale units are selected 803 as described above with reference to Fig. 2. Density- 
velocity and porosity-velocity relationships may be derived interactively 805. 
Overburden stresses are obtained from log data °07 and calibrations using crossplots are 
derived from the velocity, pressure and overburden stress 809: the Eaton trend line 
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method or the Bowers method without unloading, i.e., U = 1, may be used. The derived 
calibrations are used to predict the pore pressure, density, and/or porosity or other 
selected parameters 811. The predictions of parameters of interest may be exported 819 
or displayed on a seismic section 817. Basemaps may be produced 815 using location 
5 transforms 813. 

[0107] The flow chart of Fig. 23 shows input data 901 comprising logs (e.g., sonic, 
density, gamma ray, porosity), pressure measurements in a well (RFT, MDT, mud 
weight, LOT), velocity functions, seismic horizons, elevations (for land data) and 
bathymetry (for marine data) and seismic data. Plots of logs are obtained, selected by 
lithology or by boundaries 903. Density-velocity and porosity-velocity trends may be 
obtained 905. An overburden calculation is made 907 using log data. Seismic velocities 
are obtained between defined horizons 908 and when combined with the overburden data 
907 can be used to derive calibrations incorporating boundary constraints as well as the 
effects of unloading 909. The calibrations 909 are then used to predict attributes 911 that 
may then be exported 919 or displayed on a seismic section 917. As in figs. 22-23, 
location transforms may be applied 913 and basemaps produced. 

[0108] The foregoing description has been limited to specific embodiments of this 
20 invention. It will be apparent, however, that variations and modifications may be made 

to the disclosed embodiments, with the attainment of some or all of the advantages of the 
invention. Therefore, it is the object of the appended claims to cover all such variations 
and modifications as come within the true spirit and scope of the invention. 
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